### Loading required packages ###
library(haven)
library(dplyr)

########################
##### Reading data #####
########################

## Read data - below insert the path on your computer, where you have saved the file "BES_WL_plotdata_v2.dta" 
plotdata <- read_dta("C:/.../BES_WL_plotdata_v2.dta")

####################################
######### Organizing data ##########
####################################

### The following piece of code is used to generate both Figure 1 in the main text  ###
coef.matrix.satdem_did <- matrix(c(-.0014, .013, .021, .041, -.05, -.025, .186, -.097, .063, -.047, -.11, -.10, .009, -.019, .20), nr=1)

se.matrix.satdem_did <- matrix(c(.007, .007, .007, .006, .006, .005, .005, .006, .006, .006, .007, .007, .008, .008, .008), nr=1) 

### Save to data frame (satisfaction with democracy) ###
d <- data.frame(estimate_satdem_did=as.vector(coef.matrix.satdem_did), std.error_satdem_did=as.vector(se.matrix.satdem_did))

d$x <- matrix(c(1,2,3,5,6,7,8,9,10,12,14,15,16,17,18),nr=15)

### Point estimates and standard errors for trust in MPs ###\
coef.matrix.trustmp_did <- matrix(c(-.008, .016, .005, .022, -.033, -.041, .078, -.006, .029, -.07, -.06, .007, -.038, .031, .10), nr=1) 

se.matrix.trustmp_did <- matrix(c(.006, .006, .005, .006, .005, .009, .009, .005, .008, .011, .012, .012, .01, .006, .006), nr=1) 

### Save to data frame (trust MPs) ### 
d1 <- data.frame(estimate_trustmp_did=as.vector(coef.matrix.trustmp_did), std.error_trustmp_did=as.vector(se.matrix.trustmp_did))

d1$x <- matrix(c(1,2,3,5,6,7,8,9,10,13,14,15,16,17,18), nr=15)

##########################################################
##### FIGURE 1: Separate trends for Remain and Leave #####
##########################################################

### Run the next four plots together in order to reproduce Figure 1 (uncomment the "pdf" and "dev.off" lines) ###

pdf("Figure1.pdf") 

layout(matrix(c(1, 1, 2, 2,
		    3, 3, 4, 4), nrow=2, byrow=T),
heights=c(1.5,1))


#################################################
##### Plot over Satisfaction with democracy ##### 
##### for both remain and leave #################
#################################################

plot(plotdata$wave, plotdata$mean_satdem_w, type="b", col="black", xlim=c(1,19), ylim=c(0,.6), yaxt='n', xlab="", 
ylab="", main="", xaxt='n')
axis(side=1, at=c(1,2,3,4,6,7,8,9,10,11,13,15,16,17,19), labels=c("Feb/March 2014 (W1)", "May/June 2014 (W2)", "Sep/Oct 2014 (W3)", "March 2015 (W4)", 
"May 2015 (W6)", "April/May 2016 (W7)", "May/June 2016 (W8)", "June/July 2016 (W9)", "Nov/Dec 2016 (W10)", "April/May 2017 (W11)", "June 2017 (W13)", "March 2019 (W15)", "May 2019 (W16)", "Nov 2019 (W17)", "Dec 2019 (W19)"),
cex.axis=0.7, las=2)
axis(2,seq(0,.6,.1), cex.axis=.7, las=2)
minor.tick(ny=2, tick.ratio=.5, nx=1)
abline(v=10.5, lty=2, lwd=1)#Article 50 (eff. process)
text(x=10, y=.15, labels="(1)", cex=.7)
abline(v=14.5, lty=2, lwd=1)#2 Brexit plan defeats and 1 extenstion of the Article 50 period (ineff. process)
text(x=14, y=.15, labels="(2)", cex=.7)
abline(v=15.5, lty=2, lwd=1)#1 Brexit plan defeat and 1 extenstion of the Article 50 period (ineff. process)
text(x=16, y=.15, labels="(3)", cex=.7)
lines(plotdata$wave, plotdata$mean_satdem_l,type="b", col="grey")
legend('topleft',pch=21,pt.bg=c("black","grey"),legend=c("Leave","Remain"),bty='n',cex=.8)
mtext(3,text="Levels", cex=.8)
mtext(2,text="Satisfaction with democracy", line=2.5, cex=.8)

####################################################
##### Plot over Trust in Members of Parliament #####
##### for both remain and leave ####################
####################################################

plot(plotdata$wave, plotdata$mean_trustmp_w, type="b", col="black", xlab="",
ylab="", main="", xlim=c(1,19), ylim=c(0,.6), xaxt='n', yaxt='n',)
axis(side=1, at=c(1,2,3,4,6,7,8,9,10,12,15,16,17,18,19), labels=c("Feb/March 2014 (W1)", "May/June 2014 (W2)", "Sep/Oct 2014 (W3)", "March 2015 (W4)",
"May 2015 (W6)", "April/May 2016 (W7)", "May/June 2016 (W8)","June/July 2016 (W9)", "Nov/Dec 2016 (W10)","June 2017 (W12)", "March 2019 (W15)", "May 2019 (W16)","Nov 2019 (W17)", "Dec 2019 (W18)", "Dec 2019 (W19)"),
cex.axis=0.7, las=2)
axis(2,seq(0,.6,.1), cex.axis=.7, las=2)
minor.tick(ny=2, tick.ratio=.5, nx=1)
abline(v=10.5, lty=2, lwd=1)#Article 50 (eff. process)
text(x=10, y=.15, labels="(1)", cex=.7)
abline(v=14.5, lty=2, lwd=1)#2 Brexit plan defeats and 1 extenstion of the Article 50 period (ineff. process)
text(x=14, y=.15, labels="(2)", cex=.7)
abline(v=15.5, lty=2, lwd=1)#1 Brexit plan defeat and 1 extenstion of the Article 50 period (ineff. process)
text(x=16, y=.15, labels="(3)", cex=.7)
lines(plotdata$wave, plotdata$mean_trustmp_l, type="b", col="grey")
legend('topleft',pch=21,pt.bg=c("black","grey"),legend=c("Leave","Remain"),bty='n',cex=.8)
mtext(3,text="Levels", cex=.8)
mtext(2,text="Trust in MPs", line=2.5, cex=.8)

###########################################
##### Plotting DiD estimates ##############
##### for satisfaction with democracy #####
###########################################

plot(d$x, d$estimate_satdem_did, 
xlab="", ylab="Difference in differences",
    main="", ylim=c(-.35, .25), xaxt='n', yaxt='n')
axis(side=1, at=c(1,2,3,5,6,7,8,9,10,12,14,15,16,17,18), labels=c("W1 vs. W2", "W2 vs. W3", "W3 vs. W4", "W4 vs. W6",
"W6 vs. W7", "W7 vs. W8", "W8 vs. W9", "W9 vs. W10", "W10 vs. W11", "W11 vs. W13", "W13 vs. W15", "W13 vs. W16", "W15 vs. W16",
"W16 vs. W17", "W17 vs. W19"),
cex.axis=0.7, las=2)
axis(2,seq(-.35,.35,.1), cex.axis=.7, las=2)
minor.tick(ny=2, tick.ratio=.5, nx=1)
abline(h=0.0, lty=2, lwd=1)
mtext(3,text="DiD estimates", line=0, cex=.8)

## Plotting standard errors ##
arrows(d$x, d$estimate_satdem_did-1.96*d$std.error_satdem_did, d$x, d$estimate_satdem_did+1.96*d$std.error_satdem_did, length=0, angle=90, code=3)

## Adding text ##
text(x=11, y=.18, labels="Comparison #1", cex=.7)
text(x=11.5, y=.12, labels="DiD =.063 \n p<.001", cex=.7, font=3)

text(x=16.5, y=-.17, labels="Comparison #2", cex=.7)
text(x=17, y=-.22, labels="DiD =-.10 \n p<.001", cex=.7, font=3)

text(x=13, y=-.23, labels="Comparison #3", cex=.7)
text(x=12.5, y=-.16, labels="DiD =-.11 \n p<.001", cex=.7, font=3)

text(x=15.5, y=.18, labels="Comparison #4", cex=.7)
text(x=15.5, y=.13, labels="DiD =.009 \n p>.05", cex=.7, font=3)

arrows(x0=10, y0=.12, x1=d$x[9], y1=d$estimate_satdem_did[9]+.03, lwd=1, length=0.05)
arrows(x0=14, y0=-.2, x1=d$x[11], y1=d$estimate_satdem_did[11]-.03, lwd=1, length=0.05)
arrows(x0=16, y0=.07, x1=d$x[13], y1=d$estimate_satdem_did[13]+.03, lwd=1, length=0.05)
arrows(x0=16, y0=-.15, x1=d$x[12], y1=d$estimate_satdem_did[12]-.02, lwd=1, length=0.05)

#################################
##### Plotting DiD estimates ####
##### for trust in MPs ##########
#################################

plot(d1$x, d1$estimate_trustmp_did,
xlab="", ylab="Difference in differences",
    main="", ylim=c(-.35, .25), xaxt='n', yaxt='n')
axis(side=1, at=c(1,2,3,5,6,7,8,9,10,13,14,15,16,17,18), labels=c("W1 vs. W2", "W2 vs. W3", "W3 vs. W4", "W4 vs. W6",
"W6 vs. W7", "W7 vs. W8", "W8 vs. W9", "W9 vs. W10", "W10 vs. W12", "W12 vs. W15", "W12 vs. W16", "W15 vs. W16", "W16 vs. W17",
"W17 vs. W18", "W18 vs. W19"),
cex.axis=0.7, las=2)
axis(2,seq(-.35,.25,.1), cex.axis=.7, las=2)
minor.tick(ny=2, tick.ratio=.5, nx=1)
abline(h=0.0, lty=2, lwd=1)
mtext(3,text="DiD estimates", line=0, cex=.8)

## Plotting standard errors ##
arrows(d1$x, d1$estimate_trustmp_did-1.96*d1$std.error_trustmp_did, d1$x, d1$estimate_trustmp_did+1.96*d1$std.error_trustmp_did, length=0, angle=90, code=3)

## Adding text ##
text(x=10, y=.15, labels="Comparison #1", cex=.7)
text(x=10, y=.23, labels="DiD=.029 \n p<.001", cex=.7, font=3)

text(x=16, y=-.17, labels="Comparison #2", cex=.7)
text(x=16, y=-.22, labels="DiD=-.06 \n p<.001", cex=.7, font=3)

text(x=11, y=-.2, labels="Comparison #3", cex=.7)
text(x=11, y=-.25, labels="DiD=-.07 \n p<.001", cex=.7, font=3)

text(x=14, y=.17, labels="Comparison #4", cex=.7)
text(x=14, y=.11, labels="DiD=.007 \n p>.05", cex=.7, font=3)

arrows(x0=10, y0=.12, x1=d1$x[9], y1=d1$estimate_trustmp_did[9]+.03, lwd=1, length=0.05)
arrows(x0=13, y0=-.15, x1=d1$x[10], y1=d1$estimate_trustmp_did[10]-.03, lwd=1, length=0.05)
arrows(x0=16, y0=-.13, x1=d1$x[11], y1=d1$estimate_trustmp_did[11]-.03, lwd=1, length=0.05)
arrows(x0=14, y0=.06, x1=d1$x[12], y1=d1$estimate_trustmp_did[12]+.02, lwd=1, length=0.05)

dev.off()
